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SECTION  I 


INTRODUCTION 

Boundary  element  method  (B.E.M)  solutions  to  seven  representative  two- 
dimensional  heat  conduction  problems  are  presented.  The  steady  heat 
conduction  through  an  infinite  rectangular  prism  is  the  first  case 
discussed.  Each  of  the  next  six  cases  present  increasingly  more  complex 
governing  equations  and  boundary  conditions.  The  discussion  culminates  in  the 
examination  of  heat  conduction  through  an  infinite  circular  prism  with  heat 
radiated  from  its  surface.  The  fundamental  solution  is  derived  for  each  class 
of  governing  equations  and  its  implementation  is  discussed  for  each  specific 
case.  The  time  dependent  problems  are  successfully  solved  via  the  Laplace 
transform  method.  In  all  cases  excellent  results  are  achieved  using  very  few 
grid  points. 

This  report  represents  the  preliminary  research  in  using  the  B.E.M.  to 
study  the  non-linear  coupled  partial  differential  equations  associated  with 
plastic  wave  propagation.  Its  intent  is  to  expose  the  analytical  and 
implementational  difficulties  associated  with  application  of  boundary  element 
methods  to  transient  problems. 
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SECTION  II 

THE  BOUNDARY  ELEMENT  METHOD 


A  two-dimensional  linear  differential  operator  will  be  used  to  introduce 
the  concepts  behind  the  boundary  element  method.  There  is  no  loss  of 
generality  in  specifying  a  two-dimensional  operator  since  the  three- 
dimensional  development  is  identical. 

L(u)  =  *  (i] 

and  u  =  f  on  (2a 

|£  58  9  on  r?  (2b) 


where  r  =  is  the  bounding  curve  of  surface  B  (Figure  1).  If  the  idea 

that  there  is  a  fundamental  solution  associated  Equation  1  is  accepted,  then 
the  boundary  element  method  can  be  developed  following  Rizzo  (Reference  1),  the 
solution  to  Equation  1  can  be  written  as 


u(p)  -  -  /r2U(Q)|^  (p,Q)d*  +  /r^G(p  ,Q)~(Q)dX.  -  JB  <J>(  p  ,Q)G(p  ,Q)d  a  +  W(p)  (3) 


where 


w(p)  =  -  /r  f(Q)-g£(P*Q)d*  +  Jr  g(Q)G(p.Q)dJu 


Here,  p  Is  any  point  in  B,  Q  is  any  point  on  r,  dJL  is  a  differential  arc  length 

element  and  da  is  a  differential  surface  element.  An  integral  equation  for 

-^01 

the  unknown  u(Q)  and  dnw'  can  be  obtained  through  the  limiting  process 

lim  u(p)  where  P  Is  any  point  on  the  bounding  curve.  This  process  yields 
p-*-P 

cu(P)  =  -  /r  u(Q)*§(P,Q)dJl+  /riG(P,Q)|f(Q)d*  -  /Bi»(P,Q)G(P,Q)da  +  W(P)  (4) 
with 

W(P)  =  -  /rif(Q)-|(P,Q)di  +  /r^g(Q)G(P,Q)dt. 

where  c  depends  on  the  surface  roughness  and  the  singularity  contained  in  G. 

It  is  now  possible  to  find  the  unknown  functions  u(Q)  and  — (Q).  The 
substitution  of  these  functions  into  Equation  3  gives  the  solution  to 
L(u)  *  $  in  B.  Therefore,  only  knowledge  of  the  function  on  the  boundary  is 
necessary  to  find  the  functional  value  at  any  interior  point. 

The  numerical  solution  to  Equation  4,  and  subsequently  Equation  3, 
can  be  effected  by  first  descretizing  the  boundary  curve  into  N  segments, 
i.e.  I\  ,  i=l,2,3...N.  Next,  the  functions  u(Q)  and -^-(Q)  are  approximated  by 

ft  I 

suitable  polynomials  over  each  segment.  In  this  report  u(Q)  and  -^(Q)  are 

du 

approximated  by  zero  order  polynomials,  that  is,  u(Q)  and  are  assumed 

constant  over  each  boundary  segment.  The  development  for  higher  order  poly¬ 
nomials  is  similar  and  is  given  by  Brebbia  (Reference  2).  Based  on  the 
aforementioned  assumptions  Equation  4  can  be  written  as 
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Y(P)  -+  Ef(Q.)  /  -^(P.QjdJL  -  E  g(Qi )  J  G(P,Q.)d*. 

1=1  1  rj  11  1=1  r2 

N 

+  E  J  4>(P,Qi)G(P,Q  )dc 
1=1  8  1  11 

and,  therefore  can  be  reduced  to 


[A  -  IC]{u>  +  =  {Y} 


(6) 


where 


'Ll 

where  I  is  the  Identity  matrix,  and  N2  are  the  number  of  nodes  where 

and  u  are  unknown,  respectively.  Since  {u>  has  zero  elements  in  the  positions 

where  has  non-zero  elements  and  vice  versa.  Equation  6  can  be  reduced  to 


[N  -  IC  +  M]{x>  =  {y}. 


(7) 


with 


(x)  = 


5u 

an 


au 

an 


au 

an 


N-d 


un 


N-d+1 


*N 


Equation  7  represents  a  system  of  N  equation  for  N  unknowns.  Upon 
solving  for  x  Equation  3  can  be  used  to  find  u(p)  at  any  interior  point.  A 
four  point  Gauss-Legrendre  numerical  quadrature  was  successfully  used  to 
evaluate  the  off  diagonal  elements  of  A  and  M.  This  integration  scheme  was 
also  used  to  find  the  solution  at  interior  points. 

The  chief  difficulty  encountered  when  implementing  this  numerical  scheme 
is  integrating  over  the  singularity  contained  in  G  (diagonal  elements  of  A). 

In  an  effort  to  keep  the  boundary  element  code  as  general  as  possible,  this 
integration  is  carried  out  numerically.  The  numerical  procedure  used  follows 
the  approach  suggested  by  Hornbeck  (Reference  3).  The  interval  containing  the 
singularity  is  broken  into  the  two  intervals  shown  in  Figure  2;  in  this  figure 
the  singularity  is  located  at  X^.  A  24-point  Gauss-Legendre  numerical 
quadrature  is  used  in  the  interval  to  Xj  +  .001  and  a  four-point  Gauss- 
Legendre  quadrature  is  used  from  Xj  +  .001  to  X2.  Since  the  zeroes  of  the 
Legendre  polynomials  are  clustered  near  the  ends  of  the  interval,  subdividing 
the  Interval  of  interest  places  the  maximum  number  of  integration  points  near 


the  singularity.  This  approach  proved  very  accurate  In  evaluating  the 
diagonal  elements  of  A. 


Xj_  Xj+.OOl 

Figure  2.  Division  of  the  Integration  Interval  to  Facilitate 
Accurate  Integration  Around  a  Singularity 


SECTION  III 


SYSTEMS  GOVERNED  BY  POISSON'S  EQUATION 
1.  THE  FUNDAMENTAL  SOLUTION 

The  general  governing  equation  to  steady  heat  conduction  problems  is  the 
Poisson  equation.  The  fundamental  solution  to  this  class  of  problems  will  be 
derived  following  the  approach  of  Greenberg  (Reference  4).  As  a  result,  the 
definition  of  a  fundamental  solution  will  be  made  clear. 

The  formal  statement  for  this  class  of  problem  is 


L{ T)  =  V2T  -  *  =  0 


(8a) 


where  +  is  a  function  of  position,  and  B  is  the  general  two-dimensional  region 
depicted  in  Figure  3.  To  complete  the  problem  statement  the  mixed  boundary 
conditions  are  specified, 

C  T  +  C,  *  9  on  r  (8b) 

1  Z  on 
Y 


JS 


where,  C1  and  C2  are  constants,  and  0  is  a  function  of  position.  If  G  is 
defined  as  the  fundamental  solution  to  Equation  8a  then. 


2  2 

/g  G(V2T  -  $)do  =  /g  Id q  +  jo  Ida  -  /gG^do 

3x  ay  0 


Integrating  the  first  term  on  the  right-hand  side  yields 


/B5gdxdy  MGg  -  Tf]xXJ|dy  ♦  /,T  ^  d. 
Referring  to  Figure  3,  dysinpdA  =  n -id JL;  therefore. 


VS? dxdy  *  V6^  •  T^]  '-*x  +  v 

similarly. 


J  d«  -  /  [gE  -  T-f^-ndt  ♦  /  T®3{d0 

b  ay^  r  °Y  °Y  b  ay2 


Substituting  Equation  10  and  11  into  Equation  9,  noting 
that  V2T-P  =  0,  yields 


0  =  /  -  T~]di  +  /  G did cr  +  /  T^Gda  (1 

r  B  b 

Now,  If  V*G  »  6(r-r),  vrfiere  r  defines  the  position  of  any  point  in  B,  then 


T(rJ  *  /  CTf  -  sf]di  -  /  Gdd, 


This  result  Is  Identical  to  Equation  3  for  arbitrary  mixed  boundary 
conditions.  In  general  the  last  integral  in  Equation  12  is  of  the  form 


/  TL*(G)do 
B 


where  L*  is  the  adjoint  of  L.  In  this  case  L  is  self  adjoint,  that  is 
L  *  L*.  The  preceding  paragraphs  yield  the  conditions  on  G  and 
the  definition  of  the  fundamental  solution.  Stated  mathematically 


^  -  «(c-cJ 


Replacing  the  delta  function  by  a  delta  sequence  and  writing  the  Laplace 
operator  in  polar  coordinates,  the  conditions  on  G  can  be  written  as 


where 


II 

r  dr 


-  xQ)Z  +  (y  -  yQ)2 


If -1^(0)  is  finite,  then, 


-h'nr 


Since  the  second  term  goes  to  zero  in  the  limit,  the  fundamental  solution  to 
systems  governed  by  Poisson's  equation  (which  includes  Laplace's  equation)  is 
given  by 


G  = 


(13 


At  this  point  a  distinction  between  the  fundamental  solution  and  the 
Green's  function  should  be  made. 

9 


Here,  GF  is  the  Green's  function  to  the  governing  equations  and  boundary 
conditions;  G  is  the  fundamental  solution.  GB  is  a  function  added  to  the 
fundamental  solution  to  guarantee  that  the  Green's  function  satisfies  the 
boundary  conditions  necessary  to  make  the  first  term  in  Equation  12 
vanish.  This  term  vanishing  gives  the  direct  solution  to  Equation  8a  as 


T(rJ  -  /  GF$da 
^  B 


Greenberg  (Reference  4)  gives  a  more  detailed  discussion  on  this  distinction. 

2.  IMPLEMENTATION  OF  THE  FUNDAMENTAL  SOLUTION 

With  the  fundamental  solution  derived,  the  system  defined  by  Equation  6 
can  be  further  specified.  First  the  constant  c,  produced  by  the  limit  process 
of  Equation  4,  must  be  found.  After  substituting  the  fundamental  solution 
into  a'  condensed  form  of  Equation  4,  the  limit  looks  like 


1  im  u(p)  =  1  im  /  4  d*  -  i-  -ir(Q)  1  im  /  rdt-i-lim  /  4>(p,Q)l  nr  da  (14) 

p+P  p+P  Tr  ^  p+P  r  cn  p+P 


The  evaluation  of  the  first  term  on  the  right-hand  side  yields 


lim  /  7  di  =  lim  [  /  iu  +  /  iu] 
p+P  r  e-K 3  T-er  e  r 


where  e  is  the  radius  of  a  small  semicircular  path  circumventing  the 
singularity  on  the  boundary.  In  the  limit 


lim  /  d  Jl  =  /  d i  +  it; 

p+P  r  r  r  r 


therefore,  the  contribution  due  to  the  singularity  at  r  *  0  is  it.  The  second 
integral  in 

lim  /  InrdJL  =  lim  [  /  lnrdt+  /  lndi] 
p+P  r  £->0  T-  e 

The  second  term  on  the  left-hand  side  can  be  evaluated  as 

1  im  [elne  -  1  ]  =  0; 

£■*0 

hence,  the  contribution  due  to  the  singularity  is  zero.  Since  this  is  also 
true  for  the  surface  Integral  in  Equation  14,  its  evaluation  is  omitted. 

The  substitution  of  these  Integrals  Into  Equation  4  reveals  that  c  =  1/2. 

The  system  of  Equation  7  is  now  explicit,  that  is 


The  implementation  of  system  7  (Equation  7)  is  best  described  through 
the  use  of  Figure  4.  Each  row  in  system  7  (Equation  7)  is  the  result  of 
placing  the  source  point  (xj,  yj)  at  a  node  position  on  the  boundary,  then 
integrating  its  influence  around  r.  The  source  point  is  then  moved  from 
node  to  node  until  the  circuit  is  complete.  The  next  two  sections  give 
examples  of  the  implementation  of  this  fundamental  solution  in  the  boundary 
element  method. 


Figure  4.  General  Two-Dimensional  Discreted  Domain 

3.  STEADY  HEAT  CONDUCTION  THROUGH  AN  INFINITE  RECTANGULAR  PRISM 

The  steady  heat  conduction  In  an  infinite  rectangular  prism  is  governed 

by  the  two-dimensional  Laplace's  equation.  The  domain  in  this  example  is  a 

square  with  a  non-dimensional  length  of  6.  Figure  5  gives  the  domain,  the 

boundary  conditions  and  the  grid  employed  in  the  solution. 

Y 


Figure  5.  Discretized  Square  Domain  With  a  Twelve-Point  Grid 
(Laplace  Operator) 


Since  the  normal  derivative  is  not  well  defined  at  the  corners,  nodes 
are  not  place  in  those  locations.  Jawson,  et  al.  (Reference  5)  discusses 
rounding  the  corners.  Costable  (Reference  6)  discusses  the  use  of  Mellin 
transforms,  and  Liu  and  Sutton  (Reference  7)  discuss  a  double  corner  node 
method  to  negate  the  errors  introduced  by  corners. 

In  this  case,  and  all  of  the  following  cases,  the  material  constants 
are  unity.  Table  1  presents  numerical  results  for  eight  representative 
locations  in  the  domain. 

TABLE  1.  STEADY  TEMPERATURE  DISTRIBUTION  IN  A  SQUARE 


X  Y  B.E.M.  Analytical 


2 

2 

200.54 

200.0 

2 

3 

200.40 

200.0 

2 

4 

200.54 

200.0 

7 

5 

200.79 

200.0 

3 

2 

150.01 

150.0 

3 

3 

150.01 

150.0 

3 

4 

150.01 

150.0 

3 

5 

149.96 

150.0 

4.  STEADY  HEAT  CONDUCTION  THROUGH  AN  INFINITE  RECTANGULAR  PRISM  WITH 
A  UNIT  HEAT  SOURCE 

The  addition  of  a  unit  heat  source  results  in  a  Poisson  governing 
equation.  This  section  illustrates  the  solution  to  a  simple  steady  non- 
homogeneous  partial  differential  equation.  In  this  case,  the  surface 
temperatures  are  zero.  The  same  12-point  grid  was  also  used  in  the  solution 
to  this  problem.  Table  2  gives  numerical  results  at  the  8  representative 


TABLE  2.  STEADY  TEMPERATURE  DISTRIBUTION  IN  A  SQUARE  WITH  HEAT  SOURCE 

X  Y  B.E.M.  Analytical 


2 

2 

2.175 

2.172 

2 

3 

2.399 

2.397 

2 

4 

2.175 

2.172 

yw*v»v 


SECTION  IV 


SYSTEMS  GOVERNED  BY  THE  NON-HOMOGENEOUS  MODIFIED  HELMHOLTZ  EQUATION 
1.  THE  FUNDAMENTAL  SOLUTION 

The  remaining  five  systems,  either  directly  or  via  the  Laplace  transform, 
are  governed  by  the  non-homogeneous  modified  Helmholtz  equation.  The  general 
problem  statement  is  as  follows: 

V*T  -  k2T  =  ^  (15a) 


with  boundary  conditions 


CjT  +  C 


5T  _  9 
2  bn  9 


(15b) 


where  and  C2  are  constants  and  9  is  a  function  of  position, 
definition  of  the  fundamental  solution  given  in  Section  3,  and 
modified  Helmholtz  operator  is  self  adjoint,  gives  the  equation 
fundamental  solutions  as 


Following  the 
noting  that  the 
for  the 


/g  -  k2G  *  6(x  -  x0,  y  -  yQ) 


(16) 


If  the  Laplacian  Is  written  in  polar  coordinates,  the  stretching  p  =  kr  is 
Introduced  and  G  is  restricted  to  a  function  of  r  only.  Equation  16  reduces 
to  a  zero  order  modified  Bessel  equation. 


2  d2G  .  3G  2r  . 

p  77  +  ^  -  p  G  =  6(p  "  p0) 

ap 


(17) 


15 


Since  G  should  be  singular  at  r  *  0,  the  desired  solution  to  Equation  17  Is 
the  zero  order  modified  Bessel  function  of  the  second  kind. 

G  =  C3K0(kr) 

The  constant  C3  is  a  measure  of  the  strength  of  the  singularity  at  r  *  0.  It 
can  be  evaluated  by  substituting  G  into  Equation  17  and  integrating  over  an 
arbitrarily  small  disk  encompassing  the  singularity. 

lim  /  A  fks)d a  -  lim  J  k2K  (ke)da  =  lim  /  6(e)da 
J  e+0  B  0  Je->0  B  0  e->0  B 

By  definition,  the  right-hand  side  is  unity.  By  noting  KQ ( kr)  =  -  ln(kr) 
for  | r |  <  e,  the  second  integral  on  the  left-hand  side  evaluates  to 
-  C-jnk3[(kr)2  1  n(kr)-l/2(kr)2],  and  in  the  limit  goes  to  zero.  Since  G 
is  a  function  of  r  only,  -gp  is  the  only  non  zero  normal  derivative. 

C.,  lim  /  vA  (k e)d a  =  -  lim  /  C-k  K-(ke)dJl 
e+O  B  0  e-K)  J  1 

For  | r |  <  e,  K^kr)  =  1/kr;  therefore,  this  integral  evaluates  to  -C32n. 

Hence,  C3  =  1/2  and  the  fundamental  solution  to  the  general  modified  Helmholtz 
equation  is 

G=-^Ko(kr)  (18) 

Before  proceeding  with  the  numerical  procedure  outlined  in  Section  II,  the 
constant  c  of  Equation  4  must  be  determined.  Since  this  process  Is  identical 
to  that  In  Section  III,  the  details  are  omitted.  For  piecewise  smooth  surfaces 
and  G  ■  -1/2  *k  (kr),  c  ■  3/2.  Having  found  the  fundamental  solution  and  c. 


boundary  element  methods  can  now  be  used  in  the  solution  to  a  more  complex 
class  of  heat  conduction  problems. 

As  an  aside,  the  fundamental  solution  in  Equation  18  differs  from  that 
of  Rizzo  and  Shippy  (Reference  8)  by  a  negative  sign.  In  their  case,  c  =  1/2. 

The  reason  for  the  difference  is  in  their  approach  to  finding  the  fundamental 
solution.  They  follow  Boley  and  Weiner  (Reference  9)  in  that  Rizzo  and  Shippy 
find  the  solution  to  the  conduction  equation  when  subjected  to  a  delta 
function  heat  source.  They  then  Laplace  transform  the  conduction  equation 
to  get  the  modified  Helmholtz  equation  and  Laplace  transform  the  solution  to 
get  the  fundamental  solution.  Both  fundamental  solutions  give  identical 
temperature  distributions. 

2.  STEADY  HEAT  CONDUCTION  THROUGH  A  THIN  RECTANGULAR  PLATE 

Steady  heat  conduction  through  a  thin  plate  is  the  most  elementary  case  in 
the  Helmholtz  operator  class.  It  is  steady,  homogeneous  and  the  boundary 
conditions  are  constant. 

V^T  -  T  =  0, 

T  =  1  on  r. 

A  twelve-point  grid  is  used  in  the  solution  to  this  problem.  Figure  6  shows 
the  domain  and  boundary  conditions  and,  again,  all  material  constants  are 
unity.  The  results  are  given  in  Table  3. 
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Figure  6 


Discretized  Square  Domain  With  a  Twelve-Point  Grid 
(Helmholtz  Operator) 


TABLE  3.  STEADY  TEMPERATURE  DISTRIBUTION  IN  A  THIN  RECTANGULAR  PLATE 


3.  TRANSIENT  HEAT  CONDUCTION  THROUGH  AN  INFINITE  PRISM  OF  CONSTANT 
CROSS  SECTION 

The  problems  discussed  hereafter  fall  under  the  modified  Helmoltz  operator 
class  via  the  Laplace  transform.  In  this  section,  the  Laplace  transform 
approach  to  boundary  element  methods  is  introduced.  First,  the  heat  conduction 
through  circular  and  rectangular  prisms  (both  with  unit  temperature  on  the 
boundary)  is  discussed.  Then,  the  problem  of  constant  heat  production  in  an 
infinite  cylinder  is  discussed.  Finally,  the  problem  of  heat  conduction  from 
an  infinite  circular  prism  with  heat  radiated  from  its  surface  is  discussed. 

All  the  numerical  solutions  are  compared  to  the  analytic  solutions  given  by 
Carslaw  and  Jaeger  (Reference  10). 

Before  proceeding  with  numerical  examples,  the  corresponding  transform 
domain  system  must  be  developed.  The  general  heat  conduction  equation  is 


stated  as 


and 


(19a) 

(19b) 


where  and  C2  are  constants,  and  is  a  function  of  position  and  time. 
Following  Rizzo  and  Shippy  (Reference  8),  the  Laplace  transform  of  Equation 
19a  and  boundary  conditions  Equation  19b  yields 


and 


:  A  -  sT  =  0 


C1T  +  C2  Iff  ’  ° 


(20a) 

(20b) 


where  T  signifies  the  Laplace  transform  of  T  and  s  is  the  Laplace  transform 


parameter.  In  Equation  20,  it  is  assumed  T(0,r)  ■  0.  Notice  Equation  20a 

2 

Is  identical  to  Equation  15  with  ^  *  0  and  k  =  s/k.  The  frequency  domain 
system  can  then  be  solved  by  using  the  fundamental  solution  defined  by  equation 
18.  The  time  domain  solution  can  be  retrieved  through  the  Inverse  Laplace 
transform.  In  each  of  the  following  cases,  the  least  square  method  of  Schapery 
(Reference  11)  is  used  to  perform  the  inverse  Laplace  transform.  Rizzo  and  Shippy 
(Reference  8)  also  give  a  complete  discussion  of  this  method.  In  their  paper,  Rizz 
and  Shippy  use  conditions  on  the  frequency  response  for  large  frequency  and  on  the 
time  response  for  small  time  to  find  the  steady  state  constants.  The  accuracy 
of  this  method  seems  to  be  better  when  the  steady  state  constants  are  found 
numerical ly. 

4.  HEAT  CONDUCTION  THROUGH  INFINITE  PRISMS  OF  UNIFORM  CROSS  SECTION 
WITH  UNIT  SURFACE  TEMPERATURE 

The  solution  to  transient  heat  conduction  through  an  infinite  square  prism 
is  given.  In  the  time  domain  the  surface  temperature  is  unity,  this  transforms 
to  1/S.  The  frequency  domain  problem  statement  is 


o 

II 

i 

(21a) 

T  *  1/s  on  r. 

(21b) 

An  evenly  spaced  twenty  point  grid  is  used  in  the  solution  of  this  problem. 
Figure  7  gives  the  time  responses  found  with  Equation  20. 

Next,  the  solution  to  heat  conduction  through  an  infinite  circular  prism 
with  unit  surface  temperature  is  presented.  Equations  21  also  given  this 
problem,  in  this  case  however,  an  evenly  spaced  24-point  grid  was  used. 

Figure  8  presents  the  time  response  for  a  representative  point  in  the  domain. 
The  accuracy  for  all  other  points  is  comparable. 
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Temperature  on  the  Surface 


5.  HEAT  CONDUCTION  THROUGH  AN  INFINITE  CIRCULAR  CYLINDER  WITH  CONSTANT 
HEAT  PRODUCTION 


The  next  step  In  complexity  Is  produced  by  the  addition  of  constant  heat 
production.  Here,  the  Initial  and  surface  temperatures  are  zero;  heat  Is 
produced  at  a  unit  rate  per  unit  volume.  The  frequency  domain  governing 
equation  is 

TT  -  sT  =  -1 

with 


T  =  0  on  r  =  1 

This  example  shows  the  B.E.M.  can  accurately  give  the  solution  to  frequency 
domain,  non-homogeneous  modified  Helmholtz  equation  and  hence  also  give  the 
time  response.  The  same  domain  and  24-point  grid  is  used  in  this  solution. 
Figure  9  gives  the  time  response. 


6.  HEAT  CONDUCTION  THROUGH  AN  INFINITE  CIRCULAR  CYLINDER  WITH  RADIATION 
AT  THE  SURFACE 

This  final  example  discusses  the  heat  conduction  produced  in  a  cylinder 
when  heat  is  radiated  from  its  surface.  This  example  is  different  in  that 
functional  values  for  the  boundary  conditions  are  not  explicit,  rather  a 
relationship  between  the  temperature  and  its  normal  derivative  is  given. 


and 


7*T  -  sT  = 


£  +T.I 

3n  s 


(22a) 

(22b) 


Substituting  boundary  conditions  Equation  22b  into  integral  Equation  4  gives  a 
boundary  element  equation  in  the  form  of 


Equivalently,  the  fundamental  solution.  Its  normal  derivative  and  the  boundary 
conditions  can  be  redefined  as 
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Figure 


SECTION  V 


DISCUSSION 

As  seen  in  the  previous  examples,  the  accuracy  of  the  constant  element¬ 
boundary  element  method  is  excellent  for  both  steady  and  transient  problems. 
However,  this  approach  gives  poor  and  sometimes  spurious  results  for  complex 
geometries.  As  shown  by  Cruse  (Reference  12),  the  accuracy  of  the  B.E.M.  can 
be  dramatically  increased  by  employing  linear  elements.  For  simple  goemetries 
the  small  improvement  in  accuracy  produced  by  linear  elements  does  not  justify 
the  increase  in  computer  time. 

In  general,  for  a  given  accuracy  level  the  B.E.M.  solutions  are  much  more 
efficient  than  finite  difference  or  finite  elements  solutions  (Reference  13); 
however,  the  surface  integral  introduced  by  non-homogeneous  governing 
equations  can  greatly  increase  the  run  time  of  a  boundary  integral  solution. 

To  combat  this  problem,  it  is  important  to  take  advantage  of  symmetry  when 
ever  possible.  The  run  time  of  the  example  presented  in  subsection  IV. 5  was 
reduced  up  to  85%  when  symmetry  considerations  were  invoked. 

The  chief  purpose  of  this  report  is  to  explore  the  boundary  element 
solution  to  transient  problems.  Hence,  it  is  important  to  note  that  the 
time  dependent  solutions  consistently  breakdown  for  large  time.  This  is 
exemplified  by  the  solution  to  heat  conducting  through  an  infinite  square 
prism  for  large  time  (figure  11).  Rizzo  and  Cruse  (Reference  14)  point  out 
that  this  behavior  is  attributable  to  the  inverse  Laplace  transform.  Bellman, 
et  al .  (Reference  15)  shows  that  numerical  inverse  Laplace  transforms 
inherently  grow  unstable  for  large  time.  For  short  duration  events,  such  as 
elastic  or  plastic  wave  propagation,  the  instability  does  not  present  a 
problem  (Reference  16).  In  the  case  of  long  duration  events,  such  as  creep. 
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Figure  11.  Time  Response  of  the  Temperature  in  a  Circle  With  Unit 
Temperature  on  its  Surface  (For  Large  Time) 


a  time  dependent  fundamental  solution  can  be  derived  (Reference  17).  This 
approach  is  not  as  desirable  as  it  might  seem  since  the  B.E.M.  solution 
must  be  evaluated  at  each  time  step,  hence,  greatly  increasing  computer  time. 

The  building  block  fashion  in  which  this  report  is  presented  is  intended 
to  aid  in  the  development  of  transient  boundary  element  codes  by  providing 
increasingly  more  difficult  test  cases.  To  this  purpose,  it  is  suggested 
that  the  two-dimensional  codes  given  by  either  Crouch  and  Starfield  (Reference 
18)  or  Brebbia  et  al .  (Reference  19)  be  used  as  a  starting  point.  Both  these 
codes  can  be  expanded  to  accommodate  non-homogeneous,  time  dependent  problems. 


SECTION  VI 


CONCLUSIONS 

The  success  of  the  boundary  element  method  in  the  solution  of  two- 
dimensional  transient  heat  conduction  problems  has  been  reaffirmed.  A  method 
to  find  the  fundamental  solution  is  presented:  also  a  distribution  between  the 
fundamental  solution  and  Green's  function  is  discussed.  The  important  numer¬ 
ical  considerations  lie  in  the  numerical  integration  around  the  singularity 
of  the  fundamental  solution  and  in  the  instability  of  the  inverse  Laplace 
transform  for  large  time.  The  accuracy  and  efficiency  of  the  boundary  element 
method  in  the  solution  of  linear  transient  governing  equations  represent  a 
significant  improvement  over  other  available  numerical  schemes. 
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